Missing data

The default here is simply dropping missing data from analyses (listwise deletion), but does this make sense to do?

data|>
  select(-iso2c, -iso3c, -year)|>
  visdat::vis_miss() +
  ggtitle("Missing World Bank Data for 2020") 

lm(libdem_score ~ gini, data=data)|>
  summary()

Call:
lm(formula = libdem_score ~ gini, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.56930 -0.16300  0.07276  0.21980  0.39146 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.774470   0.169130   4.579 2.26e-05 ***
gini        -0.006482   0.004760  -1.362    0.178    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2604 on 63 degrees of freedom
  (114 observations deleted due to missingness)
Multiple R-squared:  0.0286,    Adjusted R-squared:  0.01318 
F-statistic: 1.855 on 1 and 63 DF,  p-value: 0.1781

Types of missingness

Our approach to missing data will depend on our theory of why it exists.

• Missing completely at random (MCAR)
  • no systematic differences between missing and non-missing data

Effects

Potentially problematic fixes

Listwise deletion: the default

• Simply dropping observations with missing data can work as long as the data are missing completely at random.

• Standard errors will be larger than they would be if you had all those observations, but the slopes are unbiased in a linear model

  • However, that assumption may not hold in a logit or probit model

Indicator method

Set the missing variable to 0 and include an indicator for missingness

data$gini_missing <- is.na(data$gini)
data$gini_zeroes <- ifelse(is.na(data$gini), 0, data$gini)

lm(libdem_score ~ gini_zeroes + gini_missing, data=data)|>
      summary()

Call:
lm(formula = libdem_score ~ gini_zeroes + gini_missing, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.56930 -0.19196 -0.02777  0.20530  0.52354 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       0.774470   0.153152   5.057 1.06e-06 ***
gini_zeroes      -0.006482   0.004310  -1.504  0.13438    
gini_missingTRUE -0.462005   0.154736  -2.986  0.00323 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2358 on 176 degrees of freedom
Multiple R-squared:  0.1989,    Adjusted R-squared:  0.1898 
F-statistic: 21.86 on 2 and 176 DF,  p-value: 3.328e-09

Single imputation

We could set missing values at their means, or could try using a regression model to predict the missing predictors using other data.

imputation_model<-lm(gini ~ gdp_pcap + pop_over_65 + libdem_score  , data=data)

data$gini_imputed <- predict(imputation_model, newdata=data)

lm(libdem_score ~ gini_imputed, data=data)|>
  summary()

Call:
lm(formula = libdem_score ~ gini_imputed, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.50675 -0.21558  0.02165  0.21800  0.50440 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.146892   0.164428   6.975 6.71e-11 ***
gini_imputed -0.019556   0.004322  -4.525 1.14e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2464 on 168 degrees of freedom
  (9 observations deleted due to missingness)
Multiple R-squared:  0.1086,    Adjusted R-squared:  0.1033 
F-statistic: 20.48 on 1 and 168 DF,  p-value: 1.138e-05

Better Alternatives

Extreme bounds

One approach is to assume nothing and just ask “how bad could it get?”

• For categorical variables, assume all the missing observations take on their maximum and minimum possible values
• For continuous variables, there is no real “maximum” but you could use the maximum observed value, or the mean + two or three standard deviations.
  • You could also look set something more narrow as a “plausible” alternative

Extreme bounds

min_gini<-ifelse(is.na(data$gini), min(data$gini, na.rm=T), data$gini)
max_gini<-ifelse(is.na(data$gini), max(data$gini, na.rm=T), data$gini)

mincase<-lm(libdem_score ~ min_gini, data=data)
observed<-lm(libdem_score ~ gini, data=data)
maxcase<-lm(libdem_score ~ max_gini, data=data)

Extreme bounds

• Advantage: minimal assumptions (especially for limited extreme bounds on limited DVs)

• Disadvantage: might overstate uncertainty, and would be plain wrong for MCAR

This really works best when you have a very small number of missing observations. If you can show that your results are similar even under the worst case scenario assumptions, then you’re in a pretty good spot.

MICE

MICE

library(mice)
mdata<-data|>
  filter(!is.na(libdem_score))|>
  select(-country, -iso2c, -iso3c, -year, 
         -gini_zeroes, 
         -gini_imputed,
         -gini_missing
         )|>
  mutate(population = log(population),
         gdp_pcap = log(gdp_pcap)
         )
mice_data<- mice(mdata, maxit = 35, m = 10, seed = 1,      
                 printFlag =FALSE
                 )

MICE

stripplot(mice_data, gini, pch = 19, xlab = "Imputation number")

MICE

Finally, we can fit the data and then pool the results to account for uncertainty in our imputed estimates:

fit <- with(mice_data, lm(libdem_score ~ gini))

summary(pool(fit))

In practice, we would probably want to include more information here. Especially if we have things like old data from a previous analysis, we could potentially vastly improve our predictions and limit the amount of uncertainty we had to deal with here.

MICE

• MICE only works for MCAR or MAR cases. In other words: if we can’t create a regression model that can “explain” the missingness, then it doesn’t solve anything

• No free lunch: we need at least some data about each row.

• More missing data = more uncertainty.

Selection problem: truncation of Y

set.seed(500)
N = 1000
educ_year = rpois(N, 15)
earnings  = educ_year * 10 + rnorm(N, 0, 50)
lfp = ifelse(earnings < quantile(earnings, .20), 0, 1)
lfp = factor(lfp, labels=c('not in labor force', 'in labor force'))
df<-data.frame(educ_year, earnings, lfp)
df_comb <- bind_rows(
  df |> mutate(group = "true model"),
  df |> filter(lfp == "in labor force") |> mutate(group = 'observed only')
)

true_model<-lm(earnings ~ educ_year ,data=df)
observed_model<-lm(earnings ~ educ_year ,data=df[which(lfp == "in labor force"),])

ggplot(df_comb, aes(x=educ_year, y=earnings)) + 
  geom_point(aes(fill = lfp), alpha=.5, shape=21, color='black') +

  geom_smooth(method='lm', col='orange', se=FALSE) +
  facet_wrap(~group) +
  xlab("education") +
  ylab('earnings') +
  ggtitle("hypothetical wage and education results") +
  theme_minimal() +
  scale_fill_manual(values = c("white", "black") )